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I.  INTRODUCTION 


Throughout  history,  the  application  of  materials  in  engineering  design  has 
posed  a  variety  of  problems.  In  the  nineteenth  century,  the  industrial  age  ush¬ 
ered  in  a  vast  increase  in  the  use  of  metals.  It  was  soon  discovered  that  structures 
made  out  of  such  materials  were  not  perfect.  Tragic  accidents  such  as  train  wrecks 
and  bridge  collapses  soon  brought  about  widespread  concern  over  the  design  of  such 
structures.  In  many  cases,  the  blame  was  correctly  attributed  to  a  poor  basic  de¬ 
sign.  Yet  it  was  gradually  discovered  that  metals  had  deficiencies  in  the  form  of 
pre-existing  flaws,  and  such  flaws  could  initiate  cracks  and  fractures,  thus  bringing 
about  failure  of  the  structure. 

The  discovery  of  these  flaws  brought  about  an  interest  in  studying  metallic 
materials,  for  it  was  felt  that  prevention  of  such  flaws  would  improve  structural 
performance.  Over  the  next  several  decades,  the  increase  in  the  understanding  of 
metallic  behavior,  combined  with  improved  production  methods,  brought  about  a 
marked  reduction  in  the  number  of  structural  failures. 

In  the  second  World  War,  a  renewed  interest  in  the  study  of  materials  came 
about  as  a  result  of  the  failure  of  several  Liberty  ships.  Investigations  into  these 
failures  revealed  that  flaws  and  stress  concentrations  were  responsible  for  the  brittle 
fractures.  In  the  next  several  years,  high  strength  materials  were  developed  in  the 
interest  of  weight  savings.  As  many  of  these  materials  have  low  fracture  toughness, 
it  was  discovered  th<t  they  would  fail  at  stresses  below  the  service  stress  they  were 
designed  for  in  the  presence  of  small  cracks.  This  occurence  of  low  stress  fracture  in 
high  strength  materials  has  brought  about  the  development  of  fracture  mechanics. 
(1). 


Although  the  majority  of  fracture  mechanics  has  been  developed  in  the  last 
few  decades,  its  beginnings  can  be  traced  back  to  the  research  of  A.A.Griffith  [2] 
in  1921.  Griffith  argued  that  in  the  case  of  uniaxial  tensile  loading  of  a  material 
containing  a  crack  perpendicular  to  the  load,  the  crack  would  propagate  and  bring 
about  catastrophic  failure  at  a  stress  below  its  tensile  strength.  By  analyzing  the 
region  surrounding  the  crack  tip  with  respect  to  a  global  energy  balance,  Griffith 
developed  the  concept  that  a  pre-existing  crack  can  only  extend  catastrophically 
when  the  amount  of  elastic  strain  energy  released  on  growth  of  the  crack  equals  or 
exceeds  the  surface  energy  of  the  newly  formed  crack  surfaces.  The  Griffith  equation 
for  the  strength  of  a  solid  in  plane  stress  containing  a  crack  of  length  2c  was  given 
as: 


where  E  is  Young’s  modulus,  and  7  is  the  surface  energy  [2-4], 

Although  Griffith’s  theory  works  well  for  purely  brittle  materials,  it  does  not 
accurately  describe  the  fracture  situation  in  ductile  materials.  Griffith  assumed 
that  all  of  the  work  done  during  the  fracture  process  goes  into  the  creation  of  new 
surfaces;  this  does  not  allow  for  any  dissipation  of  energy  by  plastic  deformation 
and  other  energy  dissipation  mechanisms.  As  a  result,  the  original  Griffith  criterion 
was  extended  by  Irwin  and  Orowan  [5,6]  to  include  the  case  of  ductile  fracture  of 
metals. 

Irwin  and  Orowan  noted  that  the  energy  required  for  a  crack  to  extend  in 
a  metal  is  much  greater  than  the  surface  energy  of  the  new  free  surfaces.  They 
proposed  that  a  plastic  work  term  should  be  added  to  the  Griffith  surface  energy 
7-  Furthermore,  they  argued  that  7P  is  much  greater  than  7.  Thus,  7p  replaced  7 
in  Griffith’s  original  equation,  Eq.  1,  as  follows  [5-7]: 


Irwin  later  introduced  the  terms  G  and  Gic  for  the  plastic  work  associated 
with  a  spontaneously  growing  crack  [3j.  When  G,  the  energy  release  rate  or  crack 
extension  force,  becomes  equal  to  the  critical  value,  G/c,  the  crack  extension  would 
then  take  place.  This  became  known  as  the  critical  energy  relase  rate  criterion  [1,7]. 
By  substituting  Gjc  for  the  surface  energy  term  2y,  Eq.  1  becomes: 


<7/  = 


(3) 


The  validity  of  Eq.  3  is  limited  to  only  linear  elastic  behavior  of  the  specimen 
having  a  crack  with  practically  no  crack  tip  plasticity.  However,  if  there  exists  a 
significant  amount  of  plasticity,  i.e.,  the  crack  tip  plastic  zone  is  large  compared  to 
the  crack  size,  linear  elastic  fracture  mechanics  no  longer  applies.  In  1961,  A. A. 
Wells  [8]  introduced  the  critical  crack  opening  displacement  (COD)  as  a  fracture 
criterion.  This  was  originally  developed  as  a  criterion  to  treat  those  materials  that 
exhibit  a  high  ductility.  Wells  stated  that  crack  extension  is  assumed  to  occur  when 
the  COD  exceeds  a  critical  value.  One  of  the  criterion’s  major  drawbacks  is  the  fact 
that  it  does  not  permit  the  direct  calculation  of  a  fracture  stress  [1,7].  The  COD  is 
found  from: 


4(1  - 

COD  =  — - -  ,c 

if  Eo  y  3 

where  u  is  Poisson’s  ratio,  cyB  is  the  yield  strength,  E  is  Young’s  modulus,  and  Kjc 
is  the  mode  I  plane  strain  fracture  toughness. 

In  1968,  J.R.Rice  [9]  introduced  the  application  of  a  path-independent  contour 
integral  to  analyze  elastic-plastic  crack  problems.  Similar  in  principle  to  the  critical 
energy  release  rate  criterion  discussed  above,  Rice’s  J-Integral  provides  a  fracture 
criterion  for  cases  where  plasticity  effects  are  not  negligible,  and  it  is  given  as  [7] 
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J  =  J  (wdy-(Ti^)ds^  (5) 

r 

where  W  —  j  Oij  de^j  (6) 

o 

where  the  contour  T  is  traversed  in  a  counterclockwise  direction  from  one  crack  face 
to  the  other,  W  is  the  strain  energy  density,  T,  is  the  traction  vector  directed  at  a 
point  on  the  contour  T,  u,  are  the  displacement  components,  and  s  is  a  measure  of 
arc  length  along  I\  The  crack  is  located  in  the  xy  —  plane  such  that  the  crack  lies 
parallel  to  the  local  x  —  axis. 

The  J-Integral  has  been  widely  used  in  non-linear  fracture  mechanics.  It  can  be 
extended  to  critical  values  which  will  characterize  the  crack  tip  field  at  conditions 
of  imminent  fracture  initiation.  In  the  linear  elastic  case,  the  crack  growth  occurs 
if  J  exceeds  a  critical  value  Jjc  which  is  analogous  to  Gjc.  However,  J- Integral  is 
regarded  as  a  more  general  criterion  in  that  it  is  capable  of  handling  both  elastic 
and  elastic-plastic  fracture  situations  [1,10]. 

The  independence  of  J  on  the  contour  path  T  chosen  has  been  the  subject 
of  a  great  deal  of  research.  Many  results  have  indicated  that  J  tends  to  decrease 
as  the  crack  tip  contour  shrinks  to  the  crack  tip.  In  the  elastic  case  of  true  path 
independency,  however,  J  remains  constant  as  the  contour  shrinks.  In  addition, 
J  can  be  applied  to  only  stationary  cracks.  Another  limitation  is  that  the  plastic 
region  around  the  crack  tip  must  be  small  with  respect  to  the  size  of  the  region 
in  which  J  controls  the  stress  field  (lj.  In  the  opinion  of  G.C.Sih,  “There  are 
too  many  fundamentally  unresolved  difficulties  concerning  the  assocation  of  J  with 
ductile  fracture.”  [7] 

At  roughly  the  same  time  as  Rice’s  J-Integral  was  introduced,  another  fracture 
criterion  was  developed  using  the  crack  resistance  force  R  [7].  The  criterion  was 
based  on  comparing  R  to  the  energy  release  rate  G.  It  was  stated  that: 


.V 

•  d  w 


.  "V* 
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a)  for  G  <  R,  no  crack  growth 

b)  for  G  =  R,  stable  crack  growth 

c)  for  G  >  R,  fracture  instability  occurs 

An  evaluation  of  this  fracture  criterion  is  possible  upon  analysis  of  the  R-curve,  a 
plot  of  G  and  R  versus  the  crack  length.  At  this  stage,  however,  the  theory  for  the 
R-curve  has  not  been  firmly  established  [l]. 

In  1973,  G.C.  Sih  developed  a  new  criterion  based  on  the  strain  energy  in  the 
material,  and  termed  it  the  Strain  Energy  Density  Criterion  [11,12].  The  strain 
energy  density  in  a  solid  can  be  calculated  from  either  the  area  under  the  true 
stress-true  strain  diagram,  or  from  [13] : 


dXV 

dV 


1 

2  E 


+  #  <7> 


°l  +  CTy  4-  o2x  -  2 i/[axOy  +  oyoz  +  oxaz)  \ 

where  ox,  ay ,  oz  are  the  normal  stresses,  rIy  is  the  normal  shearing  stress,  E  is 
Young’s  Modulus,  G  is  the  shear  modulus,  and  v  is  Poisson’s  ratio.  From  this 
expression,  Sih’s  strain  energy  density  factor  S  can  be  analytically  derived,  and  is 
given  as: 

s  =  '(^)  M 

The  Strain  Energy  Density  Criterion  is  given  in  two  parts: 
a.  Crack  initiation  takes  place  in  a  direction  9 o  determined  by  the  relative 
minimum  of  the  strain  energy  density  factor  S: 


S  — *  Smin  at  9  =  #0 


(9) 


b.  Rapid  crack  growth  occurs  when  the  minimum  strain  energy  density  factor 
Smin  reaches  a  critical  value: 


Smin  -  Fc  (10) 

The  general  expression  for  Sc,  the  critical  strain  energy  density  factor,  is  derivable 
from  experiments  only,  and  is  given  as: 
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_  /  lirr  x  . 

Se  =  rc(~jy)c  (U) 

where  (^r)c  is  the  experimentally  determined  value  of  ^£L  at  fracture,  and  rc  is 
some  critical  radius  of  a  core  region  surrounding  the  crack  tip.  Within  the  limits 
of  linear  elastic  fracture  mechanics,  Sc  can  be  related  to  the  mode  I  plane  strain 
fracture  toughness  Kie  as  follows: 


(1  + 1/)  (1-2*4 


The  ultimate  goal  in  our  research  was  to  find  a  suitable  fracture  mechanics 
criterion  that  could  be  extended  into  a  fracture  mechanics  theory  for  composite 
materials.  In  composites,  it  is  most  likely  that  a  crack  will  not  propagate  in  a  self- 
similar  manner.  All  fracture  mechanics  theories  except  that  of  Sih’s  assume  directly 
or  indirectly  that  the  crack  must  propagate  in  a  self-similar  manner.  This  is  the 
main  reason  why  we  were  attracted  to  Sih’s  Strain  Energy  Density  Criterion  in  our 
research. 


The  main  criticism  of  Sih’s  criterion  lies  in  the  fact  that  although  the  definition 
of  5  is  a  purely  analytical  expression,  the  definition  of  Sc  is  not  purely  analytical. 
Although  r  can  be  substantiated  analytically,  rc  cannot  be;  furthermore,  rc  has  no 
physical  significance.  Likewise,  although  ^  is  well  founded  analytically,  {~$p)c  can 
be  determined  only  from  experimental  true  stress  strain  curves. 

Finite  element  analysis  was  chosen  as  the  tool  through  which  our  research  was 
performed.  One  of  the  objectives  in  our  research  was  to  determine  how  the  strain 
energy  density  behaves  in  the  presence  of  a  crack  on  a  local  level.  This  local 
is  distinguished  from  a  global  ^  in  that  the  latter  is  determined  from  calculating 
the  area  under  an  experimental  true  stress-true  strain  diagram.  Finite  element 
analysis  was  chosen  to  conduct  a  local  energy  analysis,  as  it  is  obviously  impossible 
to  conduct  such  an  analysis  experimentally. 
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Acceptance  for  employing  the  finite  element  method  in  non-linear  situations 
is  dependent  on  several  factors.  To  begin  with,  considerable  computing  power  is 
required  to  solve  problems  with  this  degree  of  complexity.  Yet  improvements  in 
digital  computers  over  the  last  decade  of  so  have  increased  computing  power  while 
lowering  the  computing  costs.  Secondly,  the  accuracy  of  any  proposed  solution 
technique  must  be  proven  before  non-linear  finite  element  analysis  can  be  applied 
to  design. 

In  the  present  study,  an  elastic-plastic  stress  analysis  was  performed  on  A-517 
steel  plates  in  order  to  observe  the  local  behavior  of  the  strain  energy  density  in 
the  vicinity  of  the  crack.  A-517  is  a  low  carbon,  quenched  and  tempered  alloy  steel 
intended  primarily  for  use  in  welded  bridges  and  other  structures.  Compared  to 
other  high-strength  steels,  A-517  exhibits  a  very  high  yield  strength  as  well  as  good 
low  temperature  toughness  [14,15]. 

Thus,  the  objectives  of  the  present  investigation  were  as  follows: 

1.  Develop  a  suitable  finite  element  analysis  program  to  perform  a  non-linear 
stress  analysis  of  cracked  solids. 

2.  Use  this  finite  element  program  to  study  local  in  cracked  A-517  steel 
plates. 

The  finite  element  method  is  a  numerical  analysis  technique  utilized  for  ob¬ 
taining  solutions  to  an  extensive  variety  of  engineering  problems.  Any  continuous 
quantity,  such  as  stress,  strain,  pressure,  or  temperature,  can  be  approximated  by 
models  constructed  of  a  set  of  piecewise  continuous  functions  defined  over  a  finite 
number  of  elements.  Using  the  concept  that  any  continuous  function  can  be  rep¬ 
resented  by  linear  combinatio  \s  of  algebraic  polynomials,  approximation  functions 
are  derived  for  each  of  these  elements.  Elements  are  connected  at  common  nodal 
points  and  collectively  approximate  the  shape  of  the  domain.  The  assemblage  of 
elements  is  based  on  the  concept  that  the  solution  is  continuous  at  the  boundaries 


common  to  the  elements  [16-20]. 

There  are  generally  six  steps  followed  by  the  computer  in  the  finite  element 
method: 

1.  Read  in  input  data  (including  the  idealization) 

2.  Select  interpolation  or  approximation  function 

3.  Compute  the  properties  desired  for  each  element 

4.  Assemble  the  element  properties 

5.  Obtain  the  system  of  equations 

6.  Solve  the  system  of  equations 

The  first  step  in  preparing  the  input  data  is  setting  up  an  idealization  to  rep¬ 
resent  the  structure  (or,  more  generally,  the  domain).  The  resulting  solutions  for 
the  program  will  obviously  depend  on  the  idealization  created;  it  is  here  where 
experience  in  using  finite  element  analysis  is  evidenced.  Once  certain  guidelines 
are  established,  however,  even  a  beginner  can  develop  effective  idealizations.  These 
guidelines  are  as  follows:  [21] 

1.  Lay  out  the  structure  to  scale,  preferably  on  linear  graph  paper 

a.  If  the  structure  is  symmetrical,  situate  the  structure  so  that  one  of  the 
coordinate  axes  corresponds  to  the  axis  of  symmetry. 

b.  Use  enlargements  or  blowup  areas  where  necessary 

2.  Divide  the  structure  into  a  suitable  number  of  elements;  concentrate  the  ele¬ 
ments  in  areas  of  high  stress 

3.  Sketch  in  intermediate  nodes  along  the  element  edges 

4.  Number  the  nodes  using  a  suitable  technique 

5.  Number  the  elements  using  a  suitable  technique 

6.  Compute  nodal  coordinates 

A  word  on  various  types  of  elements  is  due  at  this  point.  The  long-used 
constant-stress  triangular  element  is  now  felt  to  be  both  obsolete,  inefficient,  and 


surprisingly  inaccurate  [22].  This  is  readily  evidenced  upon  comparison  to  to¬ 
day’s  higher  order  elements  such  as  the  12-node  isoparametric  quadrilateral  element 
(Quad-12  elem^it).  See  Table  1  and  Figure  1.  The  use  of  higher-order  elements  pro¬ 
duces  more  accurate  results  in  those  areas  where  the  gradient  cannot  be  accurately 
approximated  by  sets  of  constant  values.  As  opposed  to  the  constant  stress  trian¬ 
gular  elements,  the  Quad-12  element  has  a  continuously  varying  stress  field  across 
its  face.  The  displacement  varies  cubically  within  the  Quad-12  element,  as  opposed 
to  linearly  in  the  constant-stress  element.  As  a  result,  one  Quad-12  element  can 
replace  as  many  as  200  triangular  elements,  thus  reducing  data  preparation  time 
and  computer  CPU  time.  In  addition,  the  use  of  the  Quad-12  element  increases  the 
accuracy  of  the  solution. 

The  advantages  of  the  Quad- 12  element  are  thus: 

1.  Displacements  vary  cubically  over  the  element.  The  element  approximates  the 
true  displacement  function  with  a  third-degree  polynomial  degree  fit. 

2.  The  geometry  of  the  element  edges  may  vary  cubically;  thus,  curved  edges  may 
be  used  to  more  closely  approximate  the  structure. 

3.  As  the  stresses  are  given  by  appropriate  first  derivatives  of  the  displacement 
functions,  they  vary  quadratically  over  the  element. 

In  numbering  the  nodes  and  elements,  a  technique  is  used  in  order  to  maximize 
the  efficiency  of  the  computer  program.  As  a  major  part  of  the  solution  procedure 
in  finite  element  analysis  is  the  mathematical  manipulation  of  matrices,  the  compu¬ 
tation  time  is  directly  related  to  the  size  of  these  matrices.  The  set  of  equations  that 
arise  have  a  large  number  of  coefficients  which  are  zero.  Upon  analysis  of  a  typical 
system  matrix,  it  is  seen  th_t  all  non-zero  coefficients  will  fall  •  .;thin  two  imaginary 
lines  which  can  be  constructed  parallel  to  the  main  diagonal.  The  distance  from  the 
diagonal  to  this  imaginary  line  is  commonly  referred  to  as  the  bandwidth.  Reducing 
the  bandwidth  reduces  the  storage  requirements  for  the  program,  thus  reducing  the 


9 


total  computation  time. 

The  most  effective  way  to  minimize  the  bandwidth  is  to  number  the  nodes  in  the 
structural  idealization  in  such  a  fashion  as  to  keep  the  difference  between  numbers 
of  adjacent  nodes  to  a  minimum.  Similarly,  the  elements  in  the  idealization  should 
be  numbered  in  the  same  fashion  so  as  to  minimize  the  difference  between  the 
numbers  of  adjacent  elements  [19, 21].  The  following  section  deals  with  the  problem 
formulation. 


Table  1.  Tabular  comparison  of  triangular  and  Quad-12  elements 
for  a  cantilever  beam  problem 


Factor 

Triangular 

Quad- 12 

Number  of  Nodes 

400 

243 

Number  of  Elements 

686 

3 

Semibandwidth  of  Stiffness  Matrix 

20 

24 

Order  of  Stiffness  Matrix 

800 

56 

Required  Computer  Time  on  CDC  6400,  sec. 

53 

6 

Average  Displacement  Error,  percent 

-9 

fts  0 

Average  Stress  Error,  percent 

-20 

S3  0 

Manhours  to  prepare  data 

2 

0.2 

Turnaround  Time 

overnight 

5  minutes 

H.  PROBLEM  FORMULATION 


The  first  objective  in  the  project  was  to  develop  a  suitable  finite  element  anal¬ 
ysis  program  in  order  to  study  the  local  behavior  of  in  cracked  steel  plates. 
PAPST,  a  two-dimensional  elastic-plastic  finite  element  program,  was  selected  as 
a  starting  software  package.  PAPST,  an  acronymn  for  Plastic  Axisymmetric  and 
Planar  STructurcs,  has  been  developed  by  the  Navy  over  the  last  several  years  for 
the  analysis  of  elastic-plastic  crack  problems  [21-23]. 

The  features  of  PAPST  are  numerous.  The  program  incorporates  the  12-node 
quadrilateral  isoparametric  element  adapted  to  plane  stress,  plane  strain,  and  ax- 
isymmetric  conditions  of  structural  behavior.  The  displacements  and  geometry  vary 
cubically  in  the  Quad-12  element,  thus  allowing  the  accurate  modeling  of  curved 
structural  boundaries.  By  the  nature  of  the  Quad-12  element,  relatively  few  ele¬ 
ments  are  needed  to  effectively  model  most  simple  structures.  PAPST  uses  three 
basic  elements  in  all:  the  standard  Quad-12  element,  the  small  circular  core  element 
which  completely  encloses  a  crack  tip,  and  the  enriched  12-node  isoparametric  ele¬ 
ment  which  has  a  corner  node  corresponding  to  a  crack  tip  [22]. 

The  outstanding  feature  of  PAPST  is  its  non-linear  stress  and  fracture  mechan¬ 
ics  analysis.  It  has  the  ability  to  treat  cracks  in  fracture  mechanics  applications  by 
two  different  methods,  both  of  which  will  produce  a  direct  calculation  of  the  Mode  I 
and  Mode  II  stress  intensity  factors.  Rice’s  J~Integral  can  be  evaluated  on  up  to  ten 
different  paths  surrounding  the  crack  tip.  The  J-Integral  for  elastic-plastic  crack 
problems  is  given  by  Eq.  5. 

PAPST  incorporates  the  incremental  flow  theory  of  plasticity  and  the  Von 
Mises  yield  criterion  with  isotropic  or  kinematic  hardening  laws,  or  a  combination 
of  the  two.  The  incremental  plasticity  theory  is  given  as  [23]: 
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where: 


iij=  deviatoric  strain  rate  components  (e,j  =  e,y  —  ^ppSij) 

Sij—  current  deviatoric  stress  components  (s,y  =  o»y  —  | OppSij ) 

aJ  -=  deviatoric  stress  components  measured  from  center  of  current  yield 

surface  (s[-  =  s,y  -  aty) 

ctjj  =  coordinates  in  stress  space  of  center  of  yield  surface 
<re  =  von  Mises  effective  stress  =  ^§s,ys,y 

cri  =  4  f\  s'  s'  ■ 
ffyd  =  yield  stress 

In  modeling  the  material  in  PAP  ST,  an  experimentally  determined  true  stress- 
strain  curve  can  be  approximated  using  either  a  Ramberg-Osgood  power-hardening 
model  or  a  multilinear  model.  These  models  are  graphically  depicted  in  Fig.  2  [22] 
The  mathematical  form  of  the  power-hardening  model  is  given  as:  [23] 

(15) 

(16) 

where  n,  a,  and  ay  (the  yield  stress)  are  chosen  to  best  model  the  observed  behavior. 
The  mathematical  form  of  the  multilinear  model  is  given  as:  [23] 


for 

o  <  <7y, 

.11 

a 

Cy, 

e=-  +  o 

c=  -  +  —(*!  -Oy)  +  — [02- °l)  +  ■■■+— 

where  <  o  <  <7//,  ay  is  the  yield  stress,  and 
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For  this  material  model,  the  plastic  strain  rate  is  given  by: 


tplaatic  — 


OWO~e 
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and  /(<re)  = 


ow 
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PAPST  includes  several  other  features:  [21,23] 

1.  The  crack  tip  elements  include  the  plastic  mode  I  singularity  developed  by 
Hutchinson,  as  opposed  to  the  classical  linear  elastic  singularity. 

2.  The  Newton-Raphson  iterative  procedure  is  used  to  solve  the  non-linear 
equations  under  the  restrictions  of  user-specified  convergence  criteria. 

3.  PAPST  uses  4x4  numerical  integration  (Gaussian  quadrature)  to  evaluate 
the  conventional  Quad-12  element’s  tangent  stiffness  matrices. 

4.  Loading,  unloading,  and  reloading  cycles  can  be  used  to  simulate  mechan¬ 
ically  and/or  thermally  induced  residual  stresses  and  strains. 

5.  The  first  load  increment  is  automatically  calculated  to  correspond  to  the 
first  yield  in  the  specimen;  the  following  load  incrementation  is  user- 
specified. 

6.  Simulation  of  quasi-static  crack  growth. 

7.  A  treatment  of  finite  strain  effects  through  the  use  of  an  updated  La- 
grangian  formulation. 

8.  The  capability  to  compute  the  applied  tearing  modulus  T. 

9.  Large  strain  and  large  displacement  formulation. 

As  an  outgrowth  of  the  non-linear  stress  inalysis  from  PAPST,  the  strain  en¬ 
ergy  77-  was  calculated  from  Eq.  7.  It  should  be  pointed  out  that  this  equation 
is  only  valid  in  the  elastic  region  of  the  material.  PAPST  also  has  the  capability 
to  list  those  nodes  that  have  undergone  yielding.  From  this  information,  it  can  be 
determined  where  the  formula  for  ^  can  and  cannot  be  applied.  To  overcome 


this  restriction,  can  also  be  calculated  from  the  axea  under  the  analytically 
generated  true  stress-strain  diagram.  This  approach  is  valid  for  both  the  elastic 
and  plastic  regions.  The  non-linear  stress  analysis  output  from  PAPST  can  be  used 
to  analytically  draw  the  true  stress-strain  diagram;  the  axea  under  this  curve  was 
computed  and  said  to  be  equivalent  to  local  ^r. 
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Ed.  SOFTWARE  DEVELOPMENT 


A  major  portion  of  the  research  was  devoted  to  the  development  of  a  suitable 
program  to  analyze  the  behavior  of  local  in  cracked  steel  plates.  PAPST,  a  two- 
dimensional  elastic-plastic  finite  element  program  described  extensively  in  section  II, 
was  selected  as  a  starting  software  package.  As  PAPST  was  used  mainly  for  its  non¬ 
linear  stress  analysis  capabilities,  several  modifications  and  additions  were  made  to 
the  program  to  better  serve  the  needs  of  the  present  research.  The  major  addition 
was  a  subprogram  which  performed  a  complete  energy  solution  and  analysis.  The 
modified  version  of  the  original  program  is  henceforth  referred  to  as  ARLPAPST 
in  the  rest  of  this  manuscript. 

There  were  three  major  areas  in  which  PAPST  was  modified: 

1.  Data  preparation  and  pre-processing 

2.  Modifications  to  the  running  of  the  PAPST  stress  analysis  program 

3.  Several  post-processors  added  which  perform  complete  energy  analyses 
including  both  analytical  and  graphical  interpretation  of  output  from  the 
PAPST  stress  analysis  program 

The  data  input  for  the  main  PAPST  program  is  quite  extensive.  It  includes 
nodal  connectivity  and  coordinates  (both  Cartesian  and  polar),  a  variety  of  loading 
options,  graphical  and  data  output  options,  angled  and  un-angled  crack  analysis, 
many  fracture  mechanics  analysis  options  including  calculation  of  the  J-Integral 
on  up  to  ten  user-specified  paths,  and  user-specified  non-linear  material  modeling. 
To  aid  the  user  in  entering  this  vc’uminous  amount  of  data,  an  interactive  pre¬ 
processor  was  added  to  the  original  package.  This  has  significantly  reduced  the 
data  preparation  time. 

The  second  major  area  of  modification  concerned  the  convergence  criteria  used 


by  PAPST  in  the  non-linear  stress  analysis.  The  program’s  original  convergence 
scheme  was  replaced  by  a  more  intelligent  and  versatile  one.  The  new  criteria 
reads  in  initial  convergence  parameters  specified  by  the  user  in  the  data  input.  As 
the  program  runs  through  each  load  increment,  these  parameters  are  optimized 
for  the  most  efficient  convergence  possible.  In  addition,  these  continuously  variable 
convergence  parameters  are  automatically  tested  to  determine  whether  the  resulting 
criteria  is  within  acceptable  limits. 

The  main  program  was  also  modified  to  read  in  raw  data  output  options  spec¬ 
ified  by  the  user,  thereby  printing  out  only  selected  parts  of  the  output  for  each 
increment.  These  included: 

a.  Strains  and  stresses  for  each  element 

b.  Normal  strains  at  each  node 

c.  Normal  stresses  at  each  node 

d.  Location  and  value  of  highest  Von  Mises  stress 

e.  Original  (user  and  computer  generated)  nodal  coordinates 

f.  Nodal  displacements 

The  third  and  most  extensive  area  of  modifications  to  the  main  PAPST  pro¬ 
gram  was  the  post-processing.  This  post-processing,  now  an  integral  part  of  ARL- 
PAPST,  includes  a  complete  analytical  and  graphical  energy  analysis  of  the  speci¬ 
men.  The  present  capabilities  include: 

a.  Calculation  of  area  under  the  true  stress-strain  diagram  for  any  load  in¬ 
crement  and  any  user-specified  node  in  the  mesh.  This  provides  a  direct 
means  of  calculating  local  ^r,  and  is  valid  for  nodes  in  both  the  elastic 
and  plastic  regions.  A  graph  is  produced  of  the  true  stress  strain  diagram 
in  addition  to  a  history  of  the  incremental  areas  under  the  curve. 

b.  Calculation  of  local  using  Eq.  7  for  any  node  in  the  elastic  region  at 
any  load  increment. 
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c.  Plots  of  local  (formula  or  direct)  versus  distance  of  corresponding  node 
from  crack  tip. 

d.  Plots  of  local  (formula  or  direct)  versus  0,  the  angle  of  the  correspond¬ 
ing  node  from  crack  plane. 

e.  Plots  of  the  specimen  idealization  during  loading,  indicating  the  regions 
where  yielding  has  taken  place. 

Thus,  the  original  PAPST  non-linear  stress  solver  has  been  extensively  modified 
and  improved  upon.  The  new  package,  ARLPAPST,  is  a  much  more  versatile  and 
capable  program.  Its  advantages  are  summed  as  follows: 

a.  More  user-friendly. 

b.  An  extensive  list  of  raw  data  output  and  graphical  output  options  are  now 
available. 

c.  A  more  versatile  and  efficient  convergence  criterion  for  the  non-linear  stress 
analysis. 

d.  A  complete  energy  analysis  of  the  specimen,  including  two  methods  to 
calculate  local  at  a  user-specified  node  and  at  a  user-specified  load 
increment. 

e.  Several  options  available  to  graphically  and  analytically  interpret  the  en¬ 
ergy  analysis. 


IV.  FINITE  ELEMENT  MODELING 


The  behavior  of  the  local  energy  in  the  immediate  vicinity  of  the  crack  tip 
was  analyzed  using  ARLPAPST.  All  of  the  computer  analysis  was  performed  on 
a  Digital  VAX  11/782  super  minicomputer.  The  original  specimen  analyzed  is 
pictured  in  Fig.  3.  As  it  is  symmetric  about  the  crack  axis,  only  an  idealized 
upper  symmetric  half  of  the  specimen  was  modeled.  The  idealization  shown  uses 
the  12-node  isoparametric  quadrilateral  elements  discussed  previously.  Note  how 
the  elements  are  concentrated  closer  to  the  crack  tip.  A  fine  mesh  is  required  for 
an  accurate  stress  analysis  close  to  the  crack  tip,  as  this  is  where  the  non-linear 
material  behavior  and  plasticity  is  most  evident.  But  as  the  distance  from  the 
crack  tip  increases,  the  density  of  the  mesh  can  decrease.  A  detailed  mesh  far  from 
the  crack  tip  is  both  unnecessary  and  costly. 

An  experimental  true  stress-true  strain  diagram  for  A-517  steel  was  taken  from 
the  literature  [14]  and  is  shown  in  Fig.  4.  This  curve  was  modeled  in  a  bi-linear 
fashion;  the  two  points  chosen  from  the  curve  are  represented  by: 

CTi =107,400  psi  ej  =0.00358  (yield) 

02  =  108,570  psi  £2=0.006  (post-yield) 

The  other  material  properties  used  in  the  program  were  a  Young’s  Modulus  of  30 
million  psi  and  Poisson’s  ratio  u  of  0.3.  A  load  of  6000  lb.  was  used  as  a  reference 
to  which  load  increments  were  computed  for  data  input  in  each  of  the  computer 
runs. 

Several  different  specimen  geometries  were  analyzed,  and  are  tabulated  in  table 
2.  For  each  of  these  geometries,  several  runs  were  made  to  vary  the  incremental 
loading  of  the  specimen.  When  a  new  specimen  geometry  is  first  analyzed,  a  wide 


range  of  load  increments  is  usually  chosen  in  order  to  gain  a  rough  idea  of  how 


the  specimen  responds  to  the  loading.  In  subsequent  runs,  the  load  increments  are 
concentrated  over  a  certain  range  to  more  closely  study  the  fracture  behavior  of  the 
specimen.  The  choice  of  load  incrementation  must  be  chosen  in  such  a  manner  that 
the  increments  are  close  enough  to  permit  efficient  convergence  and  an  accurate 
stress  analysis.  Yet  there  is  a  tradeoff,  for  as  the  number  of  increments  is  increased, 
so  is  the  cost  of  the  individual  computer  run.  Analytical  results  are  presented  and 
discussed  in  the  next  section. 
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Figure  4.  Experimental  True  Stress  Strain  curve  for  A-517  steel  used  for  bi-linear 
modeling 
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V.  ANALYTICAL  RESULTS  AND  DISCUSSION 


The  objective  in  our  research  has  centered  on  studying  the  behavior  of  the 
local  strain  energy  density  ^  in  the  vicinity  of  the  crack  tip.  From  the  litera¬ 
ture  it  is  found  that  ^  monotonically  decays  as  the  distance  from  the  crack  tip 
increases.  From  the  non-linear  stress  analysis  output  from  ARLPAPST,  several 
energy  analyses  were  performed  in  order  to  check  this  expected  behavior. 

A  global  analysis  of  the  energy  field  around  the  crack  tip  is  given  in  Fig.  5.  This 
is  a  global  analysis  in  the  sense  that  the  data  represents  practically  all  the  nodes 
generated  in  the  finite  element  mesh.  In  general,  it  can  be  seen  that  the  energy 
steadily  decreases  with  an  increase  in  distance  from  the  crack  tip,  as  predicted  in 
the  literature.  As  Fig.  5  indicates,  however,  there  are  some  energy  values  that 
do  not  follow  this  predicted  behavior.  This  anomalous  behavior  of  an  increase  in 
energy  away  from  the  crack  tip  was  confirmed  to  be  correct. 

To  study  this  anomalous  behavior  more  closely,  the  analysis  then  narrowed  to 
those  nodes  lying  at  a  particular  angle  9  from  the  crack  axis.  As  Fig.  6  shows,  the 
curves  at  30°  and  50°  exhibit  the  expected  behavior  of  a  decrease  in  energy  with  an 
increase  in  distance  from  the  crack  tip.  However,  the  curve  at  9  =  0°,  representing 
those  nodes  lying  directly  in  front  of  the  crack,  exhibits  a  minimum  energy  point 
followed  by  an  increase  in  energy.  This  unexpected  behavior  of  an  increase  in  energy 
with  an  increase  in  distance  only  occurs  in  those  curves  corresponding  to  9  =  0°, 
irrespective  of  the  crack  lengths  and  specimen  dimensions. 

From  Fig.  3,  it  is  seen  that  the  crack  axis  is  perpendicular  to  the  direction  of 
the  loading.  Secondly,  a  homogeneous  material  (A-517  steel)  was  used  throughout 
the  research.  Thirdly,  in  all  the  specimens  analyzed,  plane  strain  conditions  were 
maintained,  thus  eliminating  the  effects  of  plane  stress.  As  a  result,  it  is  safe  to 
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assume  that  the  crack  will  propagate  in  a  self-similar  manner  at  $  =  0°.  It  was 
felt  that  there  was  a  possible  link  between  this  anomalous  behavior  of  in  the 
0  =  0°  direction  and  the  anticipated  direction  of  crack  growth  at  9  =  0°  which  is 
coincidently  in  the  same  direction. 

The  next  step  in  the  research  was  to  study  the  effects  of  load  in  further  detail 
on  the  variation  of  the  local  strain  energy  density.  This  is  graphically  depicted  in 
Fig.  7.  As  the  distance  from  the  crack  tip  is  increased,  the  energy  initially  decreases 
rather  sharply.  The  curve  then  reaches  a  minimum  point  at  roughly  a  radius  of  0.22 
inches  from  the  crack  tip.  From  this  point  to  the  right  edge  of  the  specimen  (at  0.5 
inches  from  the  crack  tip),  an  Increase  in  energy  occurs.  At  the  lower  loads,  this 
increase  is  rather  subtle.  Yet,  at  the  higher  loads,  the  energy  increases  sharply  and 
then  appears  to  level  off. 

The  minimum  energy  point  that  exists  in  each  of  the  curves  was  referred  to  as 
( ~rpr)m»»  and  occurs  at  a  location  rm,n.  As  the  load  is  increased  on  the  specimen, 
this  value  of  (^r)mtn  becomes  more  sharply  defined,  i.e.,  the  curve  develops  a  more 
sharply  defined  cusp  at  the  bottom.  In  addition,  as  the  load  is  increased  from  1200 
lb.  to  4800  lb.,  (^r)rmn  increases.  However,  for  the  curve  at  5400  lb. 

is  significantly  lower,  and  it  signifies  that  {^r)m»n  attains  a  maximum  value  at  a 
load  between  4800  lb.  and  5400  lb.  Furthermore,  at  larger  distances  from  the  crack 
tip,  it  is  seen  that  after  an  increase  in  energy,  the  curves  tend  to  reach  a  constant 
energy  value,  a  region  in  which  the  energy  gradient  is  zero.  This  constant  value  of 
energy  was  referred  to  as  the  region  of  constant  energy  starts  occuring  at 

a  location  rcon  in  front  of  the  crack  tip. 

In  order  to  deterr'ine  the  load  at  which  (jjir)min  attains  a  maximum  value, 
ARLPAPST  was  run  with  different  load  increments  chosen  so  as  to  yield  several 
energy  values  between  the  loads  of  4800  lb.  and  5400  lb.  After  plotting  curves 
similar  to  those  in  Fig.  7,  plots  were  produced  of  the  minimum  energy  versus  the 


load.  This  is  seen  in  Fig.  8.  This  figure  shows  more  clearly  how  (^)m,n  steadily 
increases  with  an  increase  in  load  up  to  a  maximum  value  around  4980  lb.,  and  then 
sharply  decreases  with  further  increase  in  the  load.  The  peak  value  of  (^)min  was 
referred  to  as  the  load  at  which  (^)„in  occured  was  referred  to  as  a 

special  critical  load,  Pcrit- 

Our  attention  then  turned  to  studying  (^r)COn-  From  Fig.  7,  it  was  seen  that 
this  parameter  was  near  constant  at  some  load.  It  was  found  that  the  value  of 
(^r)eon  was  almost  constant,  i.e.,  the  energy  gradient  was  closest  to  zero,  at  the 
critical  load  of  4980  lb.  To  isolate  this  behavior,  Fig.  9  gives  the  energy  curve  at 
this  critical  load  taken  from  a  family  of  curves  similar  to  those  in  Fig.  7.  This  curve 
clearly  exhibits  the  very  constant  nature  of  energy  far  from  the  crack  tip  (r  >  rcon), 
and  is  designated  by  (^)*on. 

In  order  to  assess  the  significance  of  and  (^pr)*on,  it  was  imperative 

to  study  the  effects  of  different  crack  lengths  on  the  behavior  of  the  energy  fields. 
The  specimen  geometry  remained  the  same  as  mentioned  in  figure  3  (width=2  in., 
height=2.4  in.,  thickness=l  in.).  In  addition  to  the  original  crack  length  of  1.537 
inches,  specimens  with  crack  lengths  of  1.0, 1.643,  and  1.7685  inches  were  analyzed. 
It  was  found  that  a  global  energy  analysis  of  these  specimens  revealed  the  same 
anomalous  behavior  in  energy  in  the  direction  9  =  0°  as  discussed  above. 

Similar  to  the  curve  in  Fig.  8,  three  different  curves  were  drawn  in  Fig.  10  to 
further  analyze  the  behavior  of  r)mm  versus  load  for  specimens  with  different 
crack  lengths.  The  first  observation  is  that  the  peak  energy  value  ( f°r 
these  curves  are  roughly  the  same  and  are  approximately  equal  to  52  lb  —  inf  in3. 
Secondly,  beyond  the  end  of  each  curve,  it  was  found  that  plastic  collapse  occured 
in  the  specimen.  It  is  safe  to  assume  that  the  ends  of  these  curves  can  be  linearly 
extrapolated  to  cross  over  the  load  axis.  At  these  points,  (^T-)mtn  =  0;  this  would 
correspond  to  a  condition  where  all  of  the  energy  stored  in  the  specimens  has  been 
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released  and,  as  such,  the  condition  corresponding  to  the  failure  of  the  specimens. 
Thus,  the  fracture  load  can  be  analytically  predicted  for  specimens  with  different 
crack  lengths.  This  fracture  load  can  be  directly  compared  with  and  verified  by 
experimental  data.  Furthermore,  from  Fig.  10  it  was  found  that  for  all  three  curves, 
the  critical  load  corresponding  to  (^)^tn  is  approximately  83%  of  the  fracture 
load.  Thus,  signifies  a  pre-failure  condition  which  may  be  characteristic 

of  the  material. 

As  regards  the  criticality  of  [$-)c0n>  Fig.  11  and  12  represent  the  strain  energy 
density  versus  radial  distance  from  the  crack  tip  for  a  wide  range  of  external  loads; 
the  crack  lengths  used  in  these  figures  are  1.0  in.  and  1.643  in.,  respectively.  It  can 
be  seen  from  Fig.  11  and  12  that  at  a  certain  critical  load  which  is  dependent  on 
the  crack  length,  (^£r)*on  remains  constant  over  the  distance  r  >  r^r,i.  However, 
it  should  be  noted  that  the  values  of  (^r)*on  in  Fig.  11  and  12  are  the  same  and 
is  equal  to  about  280  lb  —  tn/tn3,  as  can  be  clearly  seen  in  Fig.  13,  Thus,  (^r)c0„ 
is  independent  of  the  crack  length.  It  is  also  found  that  (^)*on  is  independent  of 
the  specimen  dimensions,  as  can  be  clearly  seen  in  Fig.  14.  Recalling  that  both 
( dF-) min  an^  (w’)con  31X6  exhibited  at  a  fixed  percentage  (83%)  of  the  predicted 
fracture  loads  of  various  specimens  having  different  crack  length  and  geometry,  it 
is  observed  that  both  these  energy  quantities  refer  to  some  “pre-fracture  events” 
which  are  the  charcteristic  of  the  material.  It  is  assumed  that  (^r)J^tn  and  r)lon 
correspond  to  the  local  and  global  instability,  respectively,  prior  to  the  final  frac¬ 
ture  of  the  specimens.  (Note:  the  ratio  of  (^r)j on  to  is  5.4)  The  local 

instability  can  be  further  interpreted  as  signifying  the  crack  initiation  at  r  =  r^in, 
requiring  a  relatively  smaller  level  of  energy.  Likewise,  the  global  instability  can 
be  further  interpreted  as  signifying  a  much  higher  level  of  energy  for  possible  slow 
crack  growth,  (up  to  r  =  rCTu)  which  is  expected  in  the  A-517  steel  used  in  the 
present  investigation. 
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Figure  5.  Global  analysis  of  strain  energy  density  versus  radius  (distance)  from 
crack  tip  at  a  load  of  4800  lb. 
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Figure  6.  Comparison  of  strain  energy  density  versus  radius  for  nodes  at  3  different 
angles  from  crack  tip  at  a  load  of  4200  lb. 
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Figure  7.  Comparison  of  strain  energy  density  versus  radius  for  different  load 
increments  at  0  —  0° 
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Figure  11.  Strain  energy  density  versus  radius  curves  for  specimen  with  crack 
length  of  1.0  inches 
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Figure  13.  Effects  of  crack  length  on  strain  energy  density  versus  radius  curves 
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Figure  14.  Effects  of  specimen  width  on  strain  energy  density  versus  radius  curves 


1W.WW  .W  WWW*  'Wjvrj*rA  WJW^^wuw-jLnjUMumxju^p 


VI.  CONCLUSIONS 


*>\> 


A  computer  program  ARLPAPST  has  been  developed  to  carry  out  two-dimen¬ 
sional  non-linear  finite  element  stress  analysis  in  the  neighborhood  of  a  crack  in  A- 
517  steel  compact  tension  specimens.  A  probable  direction  of  crack  propagation  and 
loads  corresponding  to  either  catastrophic  fracture  or  plastic  collapse  of  specimens 
having  different  crack  lengths  have  been  predicted.  In  addition,  two  critical  energy 
quantities,  (^■)^un  and  on  have  been  established.  These  quantities  are  in¬ 

dependent  of  the  crack  length  and  the  specimen  geometry,  and  always  correspond 
to  a  fixed  percentage  (83%)  of  the  predicted  fracture  loads  for  different  specimens. 
For  A-517  steel,  (^pr)^,,n  and  (^pr)eon  are  52  lb  —  m/m3  and  280  lb  —  in/in 3,  re¬ 
spectively.  It  is  observed  that  (7pr)^,,„  signifies  a  local  instability  leading  to  crack 
initiation  at  a  fixed  distance  ahead  of  the  crack;  and  signifies  a  global 

instability  just  before  the  slow  crack  growth,  followed  by  the  final  failure  of  the 
specimen. 
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